#%%
import scipy
from scipy.optimize import fsolve
import numpy as np
import math
import matplotlib.pyplot as plt
from matplotlib import ticker
r_f,l_d,a,b = 0.03,0.5,4.0,32.0 

#无风险响应函数
def logic_response_fuc(r,a,b):
    return math.e**(a-b*r)/(1+math.e**(a-b*r))
#最大支付意愿密度
def wilingness_pay_fuc(r,a,b):
    return b*logic_response_fuc(r,a,b)*(1-logic_response_fuc(r,a,b))
#固定违约率下的期望利润
def e_p(r,p_b,a,b,r_f,l_d):
    return logic_response_fuc(r,a,b)*((r-r_f)*(1-p_b)-(l_d+r_f)*p_b)
#最优定价
def logic_interest_optimum(p_b,r_f,l_d,a,b):
    params = [p_b, r_f,l_d,a,b]
    def func(r, params):
        p_b,r_f,l_d,a,b = params
        return r-r_f-(1.0+math.e**(a-b*r))/b-(r_f+l_d)*p_b/(1.0-p_b)
    r = fsolve(func,[0.04], args=params)
    return r[0]
#%% 
##绘制期望利润曲线分布
fig,ax = plt.subplots()
r_list = np.linspace(0,0.4,200)
p_list = np.linspace(0,0.1,5)
for i in p_list:
    p_b_value = i
    e_list = [e_p(i,p_b_value,a,b,r_f,l_d) for i in r_list]
    plt.plot(r_list,e_list,label='$p_b=$'+str(round(i, 2)))
ax.yaxis.set_major_formatter(ticker.PercentFormatter(xmax=1,decimals=0))
ax.xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1,decimals=0))
plt.grid(True)
plt.xlabel("$interest$")
plt.ylabel("$E$")
plt.legend()
plt.show()
#%%
#绘制期望利润曲面
p_b_array = np.linspace(0.001,0.03,200)
r_array = np.linspace(0.001,0.5,200)
xx,yy = np.meshgrid(p_b_array,r_array)
zz = np.zeros((200,200),dtype=np.float32)
for i in range(0,200):
    for j in range(0,200):
        zz[i,j] = e_p(yy[i,j],xx[i,j],a,b,r_f,l_d)
fig, ax = plt.subplots(subplot_kw={'projection': '3d'})

ax.plot_wireframe(xx,yy, zz,
                  color = [0.5,0.5,0.5],
                  linewidth = 0.25)

colorbar = ax.contour(xx,yy, zz,10,
             cmap = 'RdYlBu_r')

ax.contour(xx,yy, zz, levels = [0],
           colors = '#339933',
           linestyles = '-')
fig.colorbar(colorbar, ax=ax)
ax.set_proj_type('ortho')
ax.set_xlabel('$p_b$')
ax.set_ylabel('$r$')
ax.set_zlabel('$q(r,p)$')
ax.yaxis.set_major_formatter(ticker.PercentFormatter(xmax=1,decimals=0.0))
ax.xaxis.set_major_formatter(ticker.PercentFormatter(xmax=1,decimals=0.0))
ax.zaxis.set_major_formatter(ticker.PercentFormatter(xmax=1,decimals=0.0))
plt.tight_layout()
ax.set_xlim(xx.min(), xx.max())
ax.set_ylim(yy.min(), yy.max())
ax.legend(loc='upper right', frameon=False, fontsize=15)
ax.view_init(azim=-40, elev=10)
ax.grid(False)
plt.show()

# %%
